nasapower package

This package aims at making it quick and easy to automate downloading NASA POWER(Prediction of Worldwide Energy Resource) global meteorology, surface solar energy and climatology data in your R session as a tidy data frame for analysis and use in modelling or other purposes using get_power().

Installing nasapower package

You can install the packagefrom either CRAN using:

install.packages("nasapower")

or by using the devtools package directly from github using the following code:

if (!require(devtools)) {
  install.packages("devtools")
}

devtools::install_github("ropensci/nasapower")
library(nasapower)

Using get_power() to fetch data

The get_power() function has five arguments and returns a data frame with a metadata header in the current R session. It has the following arguments:

  • Community: It can be passed with either “AG”, “SB”, or “SSE”.

    AG: provides access to the agroclimatology archive, which contains industry-friendly parameters for input to crop models.

    SB: provides access to the sustainable buildings archive, which contains parameters for the building community.

    SSE: provides access to the renewable energy archive, which contains parameters very specific to assist in the design of solar and wind powered renewable energy systems.

  • temporal_average: supported values are “DAILY”, “INTERANNUAL”, “CLIMATOLOGY”.

    DAILY : the daily average of pars by day, month and year.

    INTERANNUAL: the monthly average of pars by year.

    CLIMATOLOGY: the monthly average of pars at the surface of the earth for a given month averaged for that month over the 30 year period.

  • lonlat:

    For Single point: supply a length-two numeric vector giving the decimal degree longitude an dlatitude in that order for the data to download.

    For regional coverage: supply a length-four numeric as lower left(lon,lat) and upper right(lon,lat) coordinates as lonlat = c(xmin,ymin,ymax,ymax)

    For global coverage: to get globval coverage for CLIMATOLOGY, supply “GLOBAL” while also specifying “CLIMATOLOGY” for the argument temporla_average.

  • dates: If only one date is provided, it will be treated as both the start and the end date and only a day’s values will be returned.

    When the temporal_average is set to “INTERANNUAL”, use only two year values, eg, dates = c(1983,2010).

    This argument should not be used when temporal_average is set to “CLIMATOLOGY”.

To know the different weather values from POWER provided within this function type ?get_power, and in the arguments section, click on the highlighted parameters, which goes to a page which has all the available parameters. For rainfall, we use the “PRECTOT”.

Fetching daily data for single point:

data <- get_power(community = "SSE",
          lonlat = c(134.489563,-25.734968),
          dates = c("2000-01-01","2000-05-01"),
          temporal_average = "DAILY",
          pars = "PRECTOT")
data %>% datatable(extensions = c('Scroller','FixedColumns'), options = list(
  deferRender = TRUE,
  scrollY = 350,
  scrollX = 350,
  dom = 't',
  scroller = TRUE,
  fixedColumns = list(leftColumns = 3)
))

Fetching daily data for an area:

daily_area <- get_power(community = "AG",
          lonlat = c(150.5, -28.5 , 153.5, -25.5),
          pars = "PRECTOT",
          dates = c("2004-09-19","2004-09-29"),
          temporal_average = "DAILY")
daily_area %>% datatable(extensions = c('Scroller','FixedColumns'), options = list(
  deferRender = TRUE,
  scrollY = 350,
  scrollX = 350,
  dom = 't',
  scroller = TRUE,
  fixedColumns = list(leftColumns = 3)
))

Fetching climatology data:

global data are only available for the climatology temporal_average like we discussed earlier, setting these arguments as such will fetch the global values.

climate_avg <- get_power(community = "AG",
                         pars = "PRECTOT",
                         lonlat = "GLOBAL",
                         temporal_average = "CLIMATOLOGY"
)

climate_avg %>% datatable(extensions = c('Scroller','FixedColumns'), options = list(
  deferRender = TRUE,
  scrollY = 350,
  scrollX = 350,
  dom = 't',
  scroller = TRUE,
  fixedColumns = list(leftColumns = 3)
))

Creating spatial objects from get_power()

If you require spatial objects to work with, it is simple to convert the resultant data frame from get_power() to a spatial object using the terra::rast(type = "...")

#downloading the shapefile for the world
download.file("http://thematicmapping.org/downloads/TM_WORLD_BORDERS_SIMPL-0.3.zip" , destfile="DATA/world_shape_file.zip")

#unzipping the downloaded shapefile
system("unzip DATA/world_shape_file.zip")

shape <- readOGR( 
  dsn= paste0(getwd()) , 
  layer="TM_WORLD_BORDERS_SIMPL-0.3",
  verbose=FALSE
)

tidy_shapedf <- tidy(shape, region = "NAME")
shapefile <- readOGR(dsn = paste0(getwd()),
                     layer = "ne_10m_admin_0_boundary_lines_land",
                     verbose = FALSE)
#getting the shapefiles for the world 
map <- ne_countries(returnclass = "sf")

climate_box <- split(climate_avg,climate_avg$PARAMETER)

climate_box <- lapply(climate_box, function(x){
  x["PARAMETER"] <- NULL
  x
})
climate_box <- lapply(X = climate_box, FUN= as.matrix)

#retrieving the rainfall(precipitation) data using the above made climate_box function
PRECTOT <- rast(climate_box$PRECTOT[,c(1:2,15)],
     crs = "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs",
                 type = "xyz")

#converting above raster object into a dataframe for mapping 
PRECTOT_df <- as.data.frame(PRECTOT, xy = TRUE, na.rm = TRUE)
rownames(PRECTOT_df) <- c()

#plotting the graph
ggplot() + 
  geom_raster(data = PRECTOT_df, aes(x = x, y = y, fill = ANN)) +
  geom_sf(data = map, inherit.aes = FALSE, fill = NA) + 
  scale_fill_viridis()+
  labs(title = "Rainfall in inches",
       fill = "Annual Rainfall",
       subtitle = "Annual rainfall at various parts of the world")

Using nasapower with large geographical areas

We can plot the rainfall data for large geographical areas using the shapefiles of different countries, raster cubes and finally plotting them together to show the rainfall variation throughout the region.

library("rnaturalearth")
#getting boundary values for China
CHI <- ne_states(country = "China",
                 returnclass = "sf")

#getting individual administrative boundaries 
GA <- CHI[CHI$name == "Gansu",]
XI <- CHI[CHI$name == "Xinjiang",]
#creating a raster object
r <- rast(xmin = -180,
     xmax = 180,
     ymin = -90,
     ymax = 90,
     resolution = 0.5)
#defininf the cell structure 
values(r) <- 1:ncell(r)

#retrieving only the coordinates from the GA & XI tables 
GA_coords <- crop(r,GA)
XI_coords <- crop(r,XI)

#getting the coordinated for the first administrative boundary
GA_coords <- mask(GA_coords,vect(GA))

#converting the coordinated into a dataframe with x and y values to map
GA_df <- as.data.frame(GA_coords, xy = TRUE, na.rm = TRUE)
rownames(GA_df) <- c()

#getting the china administrative boundaries shapefile
gansu_map <- geoboundaries("China","adm1")

ggplot() + 
  geom_raster(data = GA_df, aes(x = x, y = y, fill = lyr.1)) +
  geom_sf(data = gansu_map, inherit.aes = FALSE, fill = NA) + 
  scale_fill_viridis()+
  theme_minimal()+
  labs(title = "Rainfall in Gansu, China", fill = "Rainfall")

#getting the coordinated for the first administrative boundary
XI_coords <- mask(XI_coords,vect(XI))

#converting the coordinated into a dataframe with x and y values to map
XI_df <- as.data.frame(XI_coords, xy = TRUE, na.rm = TRUE)
rownames(XI_df) <- c()

#plotting them all on the map
ggplot() + 
  geom_raster(data = XI_df, aes(x = x, y = y, fill = lyr.1)) +
  geom_sf(data = gansu_map, inherit.aes = FALSE, fill = NA) + 
  scale_fill_viridis()+
  theme_minimal()+
  labs(title = "Rainfall in Xinjiang, China", fill = "Rainfall")

Let us have a look at the monthly wise data fro rainfall throughout the world. This can be done by creating a data frame which holds the names of all the months in a year, a matrix which holds all of the data for each month, and creating a for loop for conversion of each month values into a spatRaster object.

PREC <- as.matrix(subset(climate_avg, PARAMETER == "PRECTOT")[, -3])

PREC_MONTH <- vector(mode = "list", length = 12)
names(PREC_MONTH) <- colnames(PREC)[-c(1:2,13)]
#converting each month into a spatRaster object using the for loop

for (k in colnames(PREC)[-c(1:2,13)]) {
  PREC_MONTH[[k]] <- rast(PREC[, c("LON","LAT",k)],
                        crs = "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs",
                        type = "xyz")
  
}

PREC <- c(rast(PREC_MONTH))
plot(PREC)

References